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In this paper we present trispectrum estimation methods which can be applied to general non- 
separable primordial and CMB trispectra. We review the relationship between the reduced CMB 
trispectrum and the reduced primordial trispectrum. We present a general optimal estimator for the 
connected part of the trispectrum, for which we derive a quadratic term to incorporate the effects of 
inhomogeneous noise and masking. We describe a general algorithm for creating simulated maps with 
given arbitrary (and independent) power spectra, bispectra and trispectra. We propose a universal 
definition of the trispectrum parameter Tnl, so that the integrated trispectrum on the observational 
domain can be consistently compared between theoretical models. We define a shape function for the 
primordial trispectrum, together with a shape correlator and a useful parametrisation for visualizing 
the trispectrum; these methods might also be applied to the late-time trispectrum for large scale 
structure. We derive separable analytic CMB solutions in the large-angle limit for constant and local 
models. We present separable mode decompositions which can be used to describe any primordial 
or CMB trispectra on their respective wavenumber or multipole domains. By extracting coefficients 
of these separable basis functions from an observational map, we are able to present an efficient 
estimator for any given theoretical model with a nonseparable trispectrum. The estimator has 
two manifestations, comparing the theoretical and observed coefBcients at either primordial or late 
times, thus encompassing a wider range of models, such as secondary anisotropies, lensing and 
cosmic strings. We show that these mode decomposition methods are numerically tractable with 
order f operations for the CMB estimator and approximately order for the general primordial 
estimator (reducing to order in both cases for a special class of models). We also demonstrate 
how the trispectrum can be reconstructed from observational maps using these methods. 



I. INTRODUCTION 



Single field slow-roll inflationary fluctuations in the standard picture of cosmology predict a nearly scale invariant 
spectrum of adiabatic perturbations with a nearly Gaussian distribution. Hence it can be described very accurately 
by its angular power spectrum. These predictions agree well with measurements of the cosmic microwave background 
(CMB) and large scale structure, such as those provided by WMAP and SDSS. However, it remains possible that 
there exists a mechanism for generating large non-Gaussianities in the early Universe. Measurements of such non- 
Gaussianities open up the opportunity of investigating the physics of the early universe including different inflationary 
models and competing alternative scenarios. In order to study such observations, higher order correlators, beyond 
the two-point function, offer possibly the best prospects. General methods for comparing the three point correlator, 
dubbed the bispcctrum, were developed in In those papers an integrated measure of the bispectrum was 

deflned, as well as a set of formalisms for comparing, evolving and constraining the bispectrum in the case of both the 
primordial and CMB three-point correlators. In this paper we will generalise many of these methods to the four-point 
correlator which is denoted the trispectrum. We will emphasise the application of these methods to the primordial 
and CMB trispectra. The primary motivation for this paper is to develop formalisms to bring observations to bear on 
this broader class of cosmological models. We will demonstrate that despite the complexity of trispectrum estimation, 
these methods are numerically tractable given present resources, even at Planck satellite resolution. 

In order to get large non-Gaussianity we must move away from the standard single fleld slow- roll inflation 0]. 
Multifield infiation allows the possibility for superhorizon evolution. Non-Gaussianities are generated when this 
evolution is nonlinear. We can consider superhorizon behaviour as occurring in patches separated by horizons which 
evolve independently of each other. This locality in position space translates to non-locality in momentum space and 
indicates that for such models we expect the signal to peak for ^4 ^ fci, A;2, /C3. This forms the so-called local model. 
Such models have been investigated in the context of the trispectrum in [51ll3| . Since subhorizon modes oscillate and 
so average out, the only chance to have large non-Gaussianity in single field inflationary models is when all modes 
have similar wavelengths and exit at the same time. A non-standard kinetic term allows for such a possibility. Since 
the signal peaks when the modes have similar wavelengths this class of forms are known as equilateral models and 
have been investigated using the trispectrum in P^l^EU] . It should be noted that this ampliflcation of nonlinear effects 
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around the time the modes exit the horizon is not possible for slow-roll single field inflation. It has also been shown 
in [21] [22] that a large trispectrum may be generated in the ghost inflation model. These models are so-called as 
they are based on the idea of a ghost condensate, i.e. a kind of fluid with equation of state p = — p, that can flll the 
universe, and which provides an alternative method of realising de Sitter phases in the early universe. Of course there 
are other methods to generate non-Gaussianity such as having sharp features in the potential or a non-Bunch-Davies 
vacuum. Also there are models which have features that resemble the aforementioned forms in different regimes, e.g. 
quasi-single field inflation [23], or have mixed contributions, e.g. in multifleld DBI inflation [21] . 

One of the motivations for studying the four-point correlator is that it may be possible that the bispectrum is 
suppressed but still have a large trispectrum. In particular, this behaviour may be realised in quasi-single field 
inflation [53] or in the curvaton model [5S]. It also occurs in the case of cosmic strings where the bispectrum is 
suppressed by symmetry considerations [5S1 [57] . The effects of non-Gaussianity could also be detectable in a wide 
range of astrophysical measurements, such as cluster abundances and the large scale clustering of highly biased tracers. 
In |28j the possibility of using the galaxy bispectrum to constrain the local form of the trispectrum has been reviewed. 

The trispectrum, T(fci, fc2, ^3, ^4), is generally parametrised using the variable tjsil which schematically is given by 
the ratio tj^l ~ T{k,k,k,k)/P{k)^. Standard slow-roll inflation predicts tjvl ^ r/50 where r < 1 is the tensor to 
scalar ratio [11) . Such a low signal would be undetectable since it is below the level of non-Gaussian contamination that 
would be expected from secondary anisotropics tnl ~ C'(l). Using the analysis of N-point probability distribution of 
the CMB anisotropics [29], where a local non-linear perturbative model $ = $l - {^l)) + 9NL^L+0{^j^) 

is used to characterise the large scale anistropies, the constraint —5.6 x 10^ < Qnl < 6.4 x 10^ was obtained^. For the 
more general case, there is only a weak experimental bound imposed on non-Gaussianity by the trispectrum, which is 
roughly \tnl\ ^ 10* [SO]- In [SHIM] an improved constraint on t^l was presented using estimators to allow a joint flt 
of /nl and g^L using the trispectrum of WMAP5 data. However, the analysis therein included an incomplete formula 
for the CMB trispectrum due to local non-Gaussianity^. Nonetheless, the approach indicates that vast improvements 
to trispectrum constraints should be achievable in the near future. In fact, it is expected that the Planck satellite will 
be sensitive to a value of \tnl \ ^ 560 [33] . 

The analysis of the trispectrum is a computationally intensive operation. In fact only the trispectrum induced by 
the local shape has been constrained so far by CMB data. The local form is an example of a separable shape - a notion 
which we will deflne more concretely in this paper. Essentially, since the primordial trispectrum is a six dimensional 
quantity, separability means the trispectrum is the product of one dimensional functions of each of these variables. 
Exploiting this separability reduces the problem from one of C(Zmax) operations to a more manageable C'(^max)- 
special cases we get a further reduction to C(^,^^^). 

In the next section we shall describe the CMB trispectrum and its relation to the primordial equivalent. We 
will make use of a particular parametrisation of the reduced primordial trispectrum and exploit a Legendre series 
expansion in terms of one of these parameters to write an expression for the reduced CMB trispectrum which is valid 
in general. We will also outline a general correlation method for comparing different trispectra. In this section we will 
also give a formula for the kurtosis in terms of the multipoles. In section III we deflne a shape function which is a scale 
invariant form of the trispectrum. Using this function we define a shape correlator that is expected to predict closely 
the correlation between the respective trispectra. Next, we show how to decompose this shape in order to provide a 
method for visualising trispectra. We apply this visualisation to the case of the local and equilateral models which we 
describe in section IV. We also present the Sachs Wolfe limit {I < 100) for the local and constant models. In section 
V we describe how to form a mode expansion for general non-separable shapes. This provides a rigorous method to 
find a separable approximation to any shape and therefore makes analysis of the trispectrum far more tractable. This 
expansion can be performed for both the primordial and CMB trispectra. Of immediate relevance in terms of Planck 
is to find a general measure for the size of the trisectrum. This is addressed in section VI in both the primordial 
and CMB cases. It is clearly desirable to be able to reconstruct the underlying trispectrum given the data. As we 
shall describe in section VII this is a computationally intensive task, but it is tractable. We will observe here that 
there is a degeneracy in reconstruction of the primordial trispectrum, implying that only the zeroth Legendre mode 
is recoverable. Finally in section VIII we outline a method for performing CMB map simulations for given general 
bispectra and trispectra. 



^ It should be noted that for single field local inflation tJ^£ = (l/jvz,) ■ Since Jml is constrained by the bispectrum, gjvL is the quantity 

that is constrained by the trispectrum directly in this case. 
^ The formula for the reduced local CMB trispectrum has been used in place of the full local CMB trispectrum, which appears to simplify 

the analysis. 
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II. THE CMB TRISPECTRUM 
Definition of the primordial and CMB trispectra 

We are concerned with the analysis of the four-point function induced by a non-Gaussian primordial gravitational 
potential ^(k) in the CMB temperature fluctuation field. The temperature anisotropics may be represented using 
the aim coefficients of a spherical harmonic decomposition of the cosmic microwave sky, 

^(n) ^a(„y,„(n). (f) 

Lm 

The primordial potential <f> induces the multipoles aim via a convolution with the transfer functions Ai(k) through 
the relation 

— Ai{mk)Yim{k). (2) 
The connected part of the four point correlator of the aim gives us the trispectrum. In particular, 

'-^Iimil27n2l3ni^l4m4 — (^/imi ^/2^2 ^^s^'^s '^^4*^14 ) c 

= (4^)4^)5:.'. y -^^i^ ^A,,(fci)A,,(fc2)Ai3(fc3)Az,(fc4) X 

($(ki)$(k2)$(k3)$(k4))c>li™i(ki)ri,™,(k2)r/3™3(k3)rj,„,(k4), (3) 

where ki = |kj| and the subscript c is used to denote the connected component. Naively, we would define the primordial 
trispectrum as 

(<f (ki)<i>(k2)$(k3)$(k4))c - {2^lf5{k^ + ka + ks -f- k4)T4(ki, ka, kg, k4). 
Here, the four wave- vectors form a quadrilateral as shown in Figure [ij However, a more useful definition is to write 




Figure 1: Quadrilateral defined by the four wave- vectors k,;. The diagonal is represented for K. 
($(ki)$(k2)$(k3)$(k4))e = (27r)3 j d^KSik^ + k2 + K)^(k3 + k4 - K)T<J (ki , ks , kg , k4; K). (4) 

Here the delta function indicates that the diagonal K makes triangles with (ki,k2) and (k3,k4), respectively. Of 
course there are symmetries implicit in this definition of T$ - namely, that we may form triangles with different 
combinations of the vectors. In particular, 

T3,(ki,k2,k3,k4;K)= P$(ki,k2,k3,k4;K) + J d3i^'[,5(k3 - k2 - K + K')P*(ki, kg, k2, k4; K') 

+5(k4 - k2 - K + K')P<i.(ki, k4, k3, k2; K')] (5) 
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where P$ are constructed using a reduced trispectrum 7$ via 

P$(ki,k2,k3,k4;K) = 7i(ki, ka, kg, k4; K) + 7i(k2, ki, kg, k4; K) 

+ r$(ki,k2,k4,k3;K)+r*(k2,ki,k4,k3;K). (6) 

Therefore, we need only consider the reduced trispectrum T from one particular arrangement of the vectors and form 

the other contributions by permuting the symbols. 

The CMB trispectrum may also be written in a rotationally invariant way as 

The Wigner 3j symbols impose the triangle conditions on the multipole combinations {li,l2,L) and {[3,14, L). As in 
the case of the primordial trispectrum, there are implicition symmertries in this definition. In a similar manner to 
the primordial case we can write 



Thm,i,m,iMm, = Jl^'^^^' ( ri\ m2 -m) { M ) ^isuW + ('2 ^ ^3) + {h o k), 



(8) 



with 



where the factors of powers of (—1) are induced by identities of the Wigner 3j symbol. Therefore, we again need only 
consider the reduced trispectrum T from one particular arrangement of the multipoles. Indeed we need only find the 
reduced CMB trispectrum induced by the reduced primordial trispectrum. In particular, we denote 

and observe that 

+ Ti 

+ Ti 

imil4m4l3m3l2'm2 ~^ 7l4rn4li'mil3m3l2m2 * (11) 

Relation between the primordial and CMB trispectra 

In order to relate the above definitions for the primordial and CMB trispectra we use the following identities 

e"" = tej;.'i,(fc)l'i™(C)i;'„(i'), 

Yi_m = {-ITYC^. (12) 
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We find using these identities with equations ^ and Q 

Tl,m,hm2hmshrn, = H)^'' J W^fc, A;, (fc, (k, )) rf-'^i^T* (ki , , kg , k4 ; K) 

J.' J." .ml M' M" 



X 

li,L',L" m'.,M',M'' 

X [JL' (i^rOr^^M' (K)12'M' (f l)] [JL" (i^r-2)l2„M" (K)rL"M" (^2)] . (13) 



where kj represents the unit vector in the direction k,; 




Figure 2: Quadrilateral defined by the four wavenumbers ki, the diagonal K, and the angle 64 out of the plane of the first 
triangle. 

To calculate further, we must choose an appropriate parametrisation for 7$. We note that the primordial trispectrum 
shape has 6 degrees of freedom. We could define the quadrilateral uniquely by the lengths of the four sides ki = jk^l, 
together with the two diagonals K = |K| and K — |K|. However, we find it more convenient to represent the sixth 
degree of freedom with the angle 04 which represents the deviation of the quadrilateral from planarity (as illustrated 
in Figure [2]). Many well-motivated primordial models, such as the local and equilateral cases we shall discuss, are 
planar (i.e. 64^ = 0). So we choose the independent parameters to identify the shape to be (fci, ^2, ^03, ^4, iiT, ^4), that 
is, 7$ — 7$(fci, A:2) ^3; ^4! K, 94). With this parametrisation, we have the following identities. 



Att Y / \ nil 1712 

dnfYi„^{r)Y{;^,{r) = Sii,d,nm' (14) 



and ( 12 1 we find 



Ti^raihrnoMra^iira^ = ( " j ^ ^ (^l)*^ ( {kik2kzk4K)'^ dkidk2dk3,dkAdKrldrirldr2i hiKri) j L{Kr2) 

X [j,, ik,n)Ai, (fci)] [Ji. (fc2ri) A,, ik2)] [ji, {k,r2)Ai, (k,)] [j,, (^^4^2) A,, (fc4)]/i;,;.L'/i/3i^L' (-l)™^ 
^ (mi m2 -M') (i M') / r*(fci, ^2, ^3, /^4; i^, ^^4)11.™.(k4) V^^(k4)(15) 



where we write 



n2/i + l)(2?2 + l)(2£' + l) / k h L' 
~ V 47r I 



(16) 
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Next, we note that inverting equation (10) gives the expression 



mi,M 

The sum over mi,m2 is proportional to 



mi 7712 M J \ 7713 ?Ti4 — Af 



\ 777,1 Tn.T A'/ / \ 777,1 mn — A/ / ' 



TTll 77l2 Af J \ nil 7772 — Af 



and therefore the sum over L', M' implies L' = L and M' = — Af. The sum over 7773, M is then proportional to 

^ ( 777^3 777*4 -M^ { ms -m'^ -Af ) = ^j7^^i^'^4^"'i^-"^'4- (1^) 



ma, A/ 

Combining these we find that 



Ti[^i[^{L) = hi^i^Lhi,ihL (^^^ J {kik2k:ikiKYdkidk2dk:>,dkidKrldrirldr2iL{Kri)jL{Kr2) 
X [7,, (fcin) A,, (/ci)] [jj, (fc2ri) Ai, (fca)] [ji3 (A:3»-2) Aj,, (fcg)] b(4 (^47-2) A,, {k^)] 

^^z^ E j dn^T^{k,,k2MM-,KMyun^Myi:^s^i) m 

m4 — — /4 

We can decompose this expression further by expanding the primordial trispectrum as a Legendre series. In particular, 
we write 

00 

T^{ki,k2.,k3,ki;K,0i) = E "^."(^i^ ^2, ^3, ^4; Ji')-F'n(cos 6*4). (21) 

This is an expansion about the 77 = planar mode which, as we have noted, is sufficient for describing many well- 
motivated models. Noting that Pn = \J -if^i ^wO; our expression for the CMB trispectrum becomes 

Tif^{L) = hi,i,Lhi,uL j {k^k2k^k^Kfdk^dk2dk^dk^dKrldnrldr2jL{Kr^)iL{Kr2) 

X [ji^ {k,n)Ai^ (k,)] [ji^ {k2n)Ai^ (fca)] [ji, (/C3r2) A,3 (fcg)] [j^ (fc4r2) A,, (k^)] 

This expression may be further simplified by noting 



X 

m4 — — ^4 n— 



and 



000; ^ ' 

which together imply that the final line reduces to 7$,o- In particular, we have 
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rl\['{L) = hi^i^Lhi.uL 0) J{kik2k3k4Kfdkidk2dk3dk4dKrldnrldr2jL{Kn)jL{Kr2) 

X [ji^ (fciri) A,, (fci)] [ji^ (fc2ri) A,, (fc2)] [ji, {k3r2)Ai^ (ks)] [j^ (fc4r2) A,, (fc4)]r*,o(fci ,k2,k3,ki;K). (22) 
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The reduction to the n = mode in ( 22 ) shows clearly that the CMB only probes and constrains the planar component 
of the primordial trispectrum r$. In order to test theories which have general non-planar n > contributions we will 
have to use 3D data, such as 21cm surveys or large-scale galaxy distributions (as we shall discuss later). 

From equation ( 20 ) it is clear that the definition reduced trispectrum [34] includes an unnecessary geometrical 
factor hij^i^ 



2Lhi^UL and we therefore advocate the use of the true reduced trispectrum, 



hi 



(23) 



by analogy with the reduced bispectrum bi-^i^i^ 
trum. To prevent confusion, however, we refer to 



BiMjh^i^i^, where Bi^i^i^ 



represents the angle-averaged bispec- 
(L) as the 'extra'-reduced trispectrum. 



Relationship between the primordial trispectrum and other probes 



As is clear from equation (22) the CMB trispectrum depends only on the zeroth Legendre mode of the primordial 
trispectrum. Therefore, in order to break this degeneracy other probes of non-Gaussianity should be considered. As 
has been discussed in [3_5, the matter density perturbations are related to the primordial fluctuations by the Poisson 
equation via the expression 

,Sk(a) -M(fc;a)$k, (24) 
where a is the scale factor and M(fc; a) is given by 

3 fc2T(fc) 
'o 



M{k;a) = --—^D+{a), (25) 



where T(fc) is the transfer function, D^{a) is the growth factor in linear perturbation theory, fi™ is the present value 
of the dark matter density and Hq is the present value of the Hubble constant. Therefore, the primordial contribution 
to the n- point connected correlation function of matter density perturbations at a given value of the scale factor is 
given by 

(<5ki [a)5^M) ■ ■ ■ '5k„(a))c - M{k^-a)M{k2; a)... M{k,,; a)($ki3'k. • ■ ■ ^kjc (26) 

Possible probes of the matter density perturbations include galaxy surveys and the Lyman alpha forest, i.e. the 
sum of absorption lines from the Ly-a transition of the neutral hydrogen in the spectra of distant galaxies and 
quasars. There are three sources of non-Gaussianity in such surveys |36) : one primordial, one due to gravitational 
instability and the last due to nonlinear bias. 21cm observations offer another probe of non-Gaussianity which are 
less subject to the unknown galaxy bias, especially at high redshift. However, uncertainties in the neutral fraction 
replaces the uncertainties in the bias in this case. There are also complications due to redshift space distortions 
arising from peculiar velocities. Despite these drawbacks, recent advances in this area suggest that probes of the 
matter density perturbations potentially represent a powerful tool to detect non-Gaussianity and possibly break the 
degeneracy implicit in trispectrum measurements using the CMB. The study of such data involves using the full 



Legendre expansion of the primordial trispectrum as in equation (21 1. In the remainder of this paper we proceed 
to investigate the CMB trispectrum. However, many of the results presented here are straightforwardly extended to 
alternative probes of non-Gaussianity as discussed here. 



Ideal Estimator 



Unfortunately the trispectrum signal, like the bispectrum, is too weak for us to measure individual multipoles 
directly. Therefore, in order to compare theory with observations it is necessary to use an estimator that sums over 
all multipoles. Estimators can be thought of as performing a least squares fit of the trispectrum predicted by theory, 
{o-iimxO-i2m2O'i3m3.0LUmi) ci to the trispcctrum obtained from observations. The trispectrum from observations is given 
by (a°^''^^a°^''^^a°^''^^a°]',^^)c where we subtract the unconnected or Gaussian part, denoted uc, from the four point 
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function, afl, afl, afl, a?^l, 

' iimi 12^2 ^3^3 (47714 



This unconnected part is related to the observed angular power spectrum Cj 

rr°^^ — ( 1 V"l +™3 A. , A A, . A -4- f 1 Vmi +7712 /^obs ^obs 

I"iir71l"b77l2 '^'37773 '*i47n4''uC — [~ ^ ) '-^l^ Oil ,12 "tTIi ,-7772 Of., ,(4 07773 , -7774 + V ) ^ h 

^ ^^^1 ,^3 *^7ni ,—777,3 ^^2 ,^4 '^7712 , — 7774 '^Z 1 , ^4 ^7T7i , — 7774 ^^2 ,^3 ^^7772 , —7773 ^ " 

We define the estimator to be 

1 / n n n n \ ( n°^^ ^obs ^obs ^obs ^ 

C_ \"'l'ni"i27T72"/3 7773 "U7774;C ^_"ij^j„i"i2 7772 ";37773"/47774;c 



obs 



Nt 



Ci^Ci^Ci^Ci, 



by 



(27) 



(28) 



where the normalisation factor Nt is given by (see Appendix A) 



{2L + l)Ci,Ci,Ci,Ci: 



(29) 



As is clear from the earlier discussion, assuming isotropy for a given theoretical model, we need only calculate the 
reduced trispectrum, TI^^{L)^ rather than the more challenging full trispectrum {aiymiO-i2m20'i3m3ai^rni)c- 

This estimator naturally defines a correlator for testing whether two competing trispcctra could be differentiated 
by an ideal experiment. Replacing the observed trispectrum with one calculated from a competing theory we have. 



li.rrii 



(ail777i 1(27772 '^/37T73 '^'47714)0(0/^^1 '^(2 7772 '^(37773 '^/47774)c 

Ci^Ci^Ci^Ci^ 



1 V Tllll{L)T[l\[^{L) 
Nt f-l {2L + l)Ci,Ci,Ci,Cu 



(30) 



where now the normalisation Nt is defined as follows. 



Nt = 



E 



Tllll{L)Tllll{L) 



(31) 



An alternative correlator between two trispcctra, which is easier to solve numerically, is found by replacing the 
trispectra by the respective reduced trispcctra in the above definitions. Therefore, when comparing two trispcctra 
we shall use this latter definition, C{T;T')- The exact relation between the two correlators can be deduced from 
Appendix B in ref . [34] ■ 



General Estimator 



The above estimator is applicable for general trispectra in the limit where non-Gaussianity is small and the observed 
map is free of instrument noise and foreground contamination. Of course, this is an idealised case and we need to 
consider taking into account the effect of sky cuts and inhomogeneous noise. Here we follow the approach of [37] (an 
approach that is further elucidated in [38j and |39j ) . As we prove in Appendix B the appropriate form of the optimal 
estimator becomes 



^general 



/sky 

N 



^121^^2 ^Ism^ ^147714) c 



liTUi 



Iimi,l2m2 

(C a° '')l^rn3{C a° '')l^rni + 3{C )iimi, (2 7772 (C )i377l3,;4 



(32) 



where 



l^m^^l'^m'^ iP'l'-^m'-^ ^12^2 ^I'^m'^ ^I'^m'^ ci 



lim. 
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/sky is the fraction of the sky outside the mask, and where the covariance matrix C is now non-diagonal due to 
mode-mode coupUng introduced by the mask and anisotropic noise. Due to the breaking of isotropy extra terms 
have been added in order to maintain the optimahty of the estimator. The optimal estimator, in the case that the 
covariance matrix is diagonal, reads 

o _ /sky {o-limi 0-l2m2 O-hrris O'Umi ) c. [^obs obs obs obs a( -|\mi^ r r obs obs 

]\l 2^ C C C C [";inil"i2m2"i3m3"i4m4 ^) '-'il b "nil -ni2 "/gma ";4m4 

Umi '3 h 

+3{—l)"^^^"^^Si^l2S„ii-m2^hU^m3~m4CliCl^ , (33) 

where Nt is given by equation (p9| . We note also that the average of this estimator is 



p I ri n n n \ / n^^^ n*^^^ n'^^^ ^obs \ 

as expected. 

In the remainder of this paper we shall refer to the ideal estimator unless otherwise stated. However, this formula 
is important for the general implementation of the formalisms introduced here. 

Kurtosis as a measure of non-Gaussianity 

As an aside, we note that the use of non-optimal estimators may also provide useful information, e.g. as a reality 
check on these complex calculations. The kurtosis of the one point temperature distribution offers such an estimator. 
The kurtosis, (72, is defined as 

92 = — ^ '-^ - 3. (35) 

(((¥(-))')) 

As we show in Appendix C (where we also include a discussion on the skewness for completeness) the kurtosis may 
be written in the following form 

The calculation of this quantity is relatively straightforward compared to the full estimator due to the absence of 
Wigner 6j symbols in the expression. 

III. THE SHAPE OF PRIMORDIAL TRISPECTRA 
Shape function 

It is known from CMB observations that the power spectrum is nearly scale- invariant. Analysis of the bispectrum 
is performed using the shape function, which is a scale invariant form of the bispectrum. To parallel this analysis we 
wish to write a scale invariant form of the trispectrum (or in particular the trispectrum modes). Therefore, we need 



to eliminate a k scaling. Motivated by ( 22 ) we define this shape by 



STiki,k2,k3,k4,K) = (fcifc2fc3fe4) ^ x^fl{ki,k2,k3,k4;K), (37) 



where N is an appropriate normalisation factor. For clarity in what follows we note that we shall use the symbol 
Sj- when referring to the shape induced by the reduced primordial trispectrum. Of course, this choice of the shape 
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function is not unique. Another choice of shape function is 

ST{ki,k2,k3,k4,K) = ^^^^'^^3^^ — T<s,^o{ki,k2,k3,k4;K), (38) 

which has the advantage of remaining independent of the diagonal K if the underlying trispectrum has this property. 
Such a class of models are discussed further in Appendix D. Nonetheless we proceed with St as our choice of shape 
function in this paper, leaving further investigation of this issue to a future publication gl]. We should also notice 
that, since our analysis here is focused on the CMB, we have only included the zeroth mode of the Legendre expansion 



as indicated by ( 22 ). However, for more general probes of non-Gaussianity, as discussed in Section II, the full Legendre 



expansion described by equation (21 1 is required. In such a case the analysis outlined here can be applied mode-by- 
mode. Due to orthogonality of the Legendre modes, extending the study is a trivial task. 
If we rewrite the reduced CMB trispectrum in terms of the shape function, 6*7-, we have 

Tl'l'iL) = Nhi,i,Lhi,uL (^l^ J dVkSr{ki,k2,k3,ki,K)K 

xAz,(fci)Ai,(fc2)Ai3(fc3)Ai4(fc4)/i':z,,3^L(fcl,fc2,fc3,fc4,i^) (39) 

where the integral is given by 



Ihhi3ULiki,k2,k3,ki,K) = J r^r2dridr2jL{Kri)jL{Kr2)ji^{kiri)ji^{k2ri)ji^{k3r2)ju{kir2) (40) 

and dVk corresponds to the area inside the region ki,K/2 e [0, fcmax] allowed by the triangle conditions. Therefore 
the shape function is the signal that is evolved via the transfer functions to give the trispectrum today. Essentially, 
l'^ acts like a window function on all the shapes as it projects from k to ^— space, that is, it will tend to smear out 
their sharper distinguishing features. This means that the shape function ^7-, especially in the scale invariant case, 
can be thought of as the primordial counterpart of the reduced CMB trispectrum 7^' (L) before projection. 



Shape correlators 

We wish to construct a primordial shape correlator that predicts the value of the CMB correlator C(T, T'). To this 
end we should consider something of the form 

F{Sr, Sr') ^ J dVkSr{ki,k2,k3, k^, K)Sr'{ki,k2, fca, k^, X)a;(fci, ^2, fcg, k^, K) (41) 

where oj is an appropriate weight function. With this choice of weight the primordial shape correlator then takes the 
form 

^' ^ ^FiSr,Sr)F{Sr',Sr'y 

The question now is what weight function should we choose? Our goal is to choose S'^uj in /c— space such that 
it produces the same scaling as the estimator T^/((2L -|~ 1)C"*) in l~ space. Let's consider the simplest case where 
k = ki = k2 = ks = k^ = K and I — li = I2 — I3 = li = L. For primordial trispectra which are scale invariant, then 

S^ik, k, k, fc, k) (jj{k, k, k, k, k) oc w(fc, k, k, k, k). (43) 

If we work in the large angle approximation, and assume I + 1 ^ I, then we know C; cx and from the analytic 
solution for the local model which we will describe below (see equations! 59) and (63l) we have 

Tll{l)^hli^^. (44) 
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Now hill oc l^^^ ) ^^'^ ^^'^ Wigner 3 J symbol has an exact solution for which 



O; ^ yamV 3Z! (Z/2)!3 ^ VVsvrr ^^^^ 



Therefore T//(Z) cx I ^ and so 



{21 + l)Cf 



oc r^. (46) 



Hence we find that we should choose a weight function u!{k,k,k,k,k) cx The particular choice of w may 

significantly improve forecasting accuracy - by, for instance, using a phenomenological window function to incorporate 
damping due to photon diffusion or smoothing due to projection from k- to I- space, but it does not impact important 
qualitative insights. A specific choice of weight function, motivated by the choice of weight function for the bispectrum, 
is the following 

wik,,k,, fc3, K) ^ ^k, + k,+Kr{k, + k,+Kr- 



Shape Decomposition 

Given the strong observational limits on the scalar tilt we expect all shape functions to exhibit behaviour close 
to scale-invariance, so that S'j-{ki,k2,k3,ki, K) will depend only weakly on the overall magnitude of the summed 
wavenumbers. Here we choose to parametrise the magnitude of the wavenumbers with the quantity 

k=^iki + k2 + K). (48) 

k is the semi-perimeter of the triangle formed by the vectors ki, k2, K. Due to the scaling behaviour the form of the 
shape function on a cross-section is essentially independent of /c, so that 

S'r(fci, k2,k3,ki, K) = f{k)ST{ku k^k, k, K) (49) 

where ki = ki/k and K — K/k. Since we are restricted to the region where the wavenumbers (fci, A:2, K) and (fca, k^, K) 
form triangles by momentum conservation, we will reparametrise the allowed region to separate out the overall scale 
k from the behaviour on a cross-sectional slice. This four-dimensional slice is spanned by the remaining coordinates. 
Concentrating on each triangle individually we reparametrise in a similar fashion to the analysis done in [3]. For 
triangle (ki, k2, K) we have 

K = k{l-(3) 

fci = ^(1 + « + /?) 

k2 = ^(l-a + /3) (50) 



while for triangle (k^, k^, K) we have 



K = ek{l-6) 

h = y(l + 7 + '5) 

fc4 = y(l-7 + '5) (51) 



12 



where e parametrises the ratio of the perimeters of the two triangles, i.e. e = ^-^^^^^j^ . We consider 1 < e < oo. The 
different expressions for K imply that 

I- P = e{l - 5). (52) 

The conditions for triangle (fci, fc2, K) that < ki, k2, K < k imply that < /3 < 1 and — (1 — /?) < a < 1 — (3, while 
the conditions for triangle {k^^kj^^K) that < ki,k2, K < ek, along with the relationship between S and (3 and the 
requirement that e > 1, imply that —{1 — P)/e < j < (1 — /3)/e- In summary, we have the following domains, 

0<k<oo, l<e<oo, 0</3<l, -(1 -/3) < a < 1 ^ < 7 < ^— ^. (53) 

e e 

With this parametrisation we can re-write the shape function and the volume element respectively as 

S-j-{ki, k2, fca, fc4, K) — f{k)S'j-{a, fi, 7, e), dV^ — dkidk2dk^dkjidK — ek'^dkdadfSd'fde. (54) 

In order to represent the shape function graphically we can choose fixed values of e in which case the shape S becomes 
three dimensional. The particular three dimensional domain is shown in Figure |3] From the image we see how the 
particular triangles created by the wavenumbers generate the three dimensional slice for each e. We can envisage the 
four-dimensional shape by imagining an orthogonal direction for e out of the page, along which are located increasingly 
squeezed rectangular pyramids. 




Figure 3: Three-dimensional shape function domain for a fixed value of e, i.e. for a particular ratio of the perimeters of the 
triangles created by the wavenumbers {k\,k2, K) and {ks,k4,K) respectively. Note that the triangle conditions on these two 
wavenumber sets restrict them to the two tetrahedral domains illustrated (right), slices through which are mapped as shown 
into the full domain, a rectangular pyramid (left). 



IV. SEPARABLE SHAPES 
Examples: Local, equilateral and constant models 

The local model is given by the reduced primordial trispectrum 

(ki , k2 , ks , k4 ; K) = ri°i (ki , k2 , kg , k4 ; K) + (ki , ka , kg , k4 ; K) (55) 
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where 

7^°i(ki,k2,k3,k4;K) = -TNLP^{K)P^{k,)P^{k3) (56) 
ri°i(ki,k2,k3,k4;K) = gNL[Pi{k2)Pi{k3)P^{ki) + Pi{ki)Piik2)Pi{ki)] . (57) 

For single field inflation we have t^l = {^Inl}'^- This relationship breaks down for multifield inflation (see 02] )■ 
We can see clearly here that the local trispectrum is independent of the angle O4, i.e. the zeroth mode of the local 
trispectrum is exactly the full local trispectrum. The primordial shapes for each of these expressions may be shown 
visually using the prescription described in the previous section and they are shown in Figures |4] and [5j As expected 
for the local model the signal peaks in the corners. However, as is easily observable the 'peaking' behaviour is 
somewhat orthogonal between the two models. Working in the Sachs- Wolfe approximation, where we replace the 
transfer function with a Bessel function, 



A/(fc) = ^ji((To - Tdoc)fc), 

the integral for the reduced trispectrum can be expressed in closed form. Setting P^{k) = A^k~^ we find 



'hl4,A 



where 



and we have used 



Similarly, 



25TNL A| ' 
9 34 



9 34 



hi^i^Lhi^UL j rldrirldr2K ^dXfc^ ^dfci/cg ^d/cg//^ (2, ri)//j2, ri) 
hi,i,Lhl,UL[\) /L(-l,l)/i,(-l,l)//3(-l,l) 



25TjVL 2^ a 
36 34 U 



2L{L + l)2h{h + l)2h{h + I) 



Ii{p,x)^ I kPdkji{k)ji{xk), 



2r^ 

1 

21(1 + !)■ 



n-hhloc/T^ _ ^ TT^ A| /2 



1 



2^2(^2 + 1)2^3(^3 + 1)2^4(^4 + 1) 



+ (^1 ^ h) 



(58) 



(59) 

(60) 

(61) 
(62) 

(63) 



Next we propose a constant model for the primordial trispectrum, analogous to the simplest model for the bispec- 
trum. This is given by 



1 



A|7V 



{kik2k3k4f K%.fiiki, k2, ks, k4; K) = 5r(fci, ^2, ^3, ^4, ^) = 1 



(64) 



with N the normalisation factor of equation (64) and the choice A|, motivated by comparison to the local model. 
Again, as for the local model, the primordial trispectrum is already a zero mode quantity with respect to angle 64, 
i.e. T = To- Using the Sachs Wolfe approximation the integral (22) can be written as 



7-/1 ^2C011St / 7- \ 



A^iV /2 



34 



-) hi^i^Lhi^UL j dxx drirJi^{Q,ri)Ii^{Q,ri)I^{Q,rix)Ii^{Q,rix)lL{l,x) (65) 




Figure 4: Local A model 1 59 1. The peak as /5 — ^ 1 corresponds to — > 0, i.e. the 'doubly-squeezed' limit. The other peak 
corresponds to fci — > 0, fca — >■ (e — l)k. As e rises above unity (i.e. for triangle {ki,k4,K) bigger than triangle (k\,k2, K)) we 
expect this peak to become suppressed as observed. 




Figure 5: Local B model (63 1. The peaks at a = 1 correspond to k2 — >■ 0, while the peaks at a = — 1 correspond to fci — > 0. For 
e = 1 (i.e. equal triangle sizes (fci, k2, K) and (^3, k^, K)) we see a more confined peaking at 7 = —1, 7 = 1, i.e. ^3 — )■ 0, ^4 — >■ 
respectively. We observe the peaking of the local B models to be somewhat orthogonal to the local A model. 
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where we write r2/ri = x. Now we can evaluate 



IiiO,x) 



2(2/ + 1)' 

TT 

2{2i + iy 



for X > 1 
for X < 1 



(66) 



2r(L + 3/2) 4' ' 2' 

7rr(L + l) r ,1 3 
2f(W)^^^^^2'^ + ^=^+2 = -) 



for a; > 1 
for a; < 1 



(67) 



where 2^1 is a generahsed hypergeometric function. We can write 2F1 in terms of a series expansion in the form 

(a)„(&)„ z" 



2Fi(a,6;c;2;) = ^ 



n=0 



(c)„ n! 



(68) 



where (p)„ = r(p + n)/T(j)). The conditions for convergence, namely that this series converges for c a non-negative 
integer with \z\ < 1, are satisfied in this case. Using this decomposition we find 



v. 



7 1/2 const 



^3/4 



1 



34 



E 



r(l/2 + n)r(L + l + n) 1 



1 

n (2/1 + \){2h + 1)(2;3 + l)(2/4 + 1) ^0 r(L + 3/2 + n) 



2n + 3 + /3 + /4 + i AiJ 2n + l + li + l2 + L \J2k ^1 



2n + 3 + Zi +^2 +i VE^« +4 ^2/ 2n + l + l3 + U + L yj^h ^ 



(69) 



Notice that this sum is still finite if the denominators Ai = li + I2 ~ I3 — I4 + 2 or A2 — I3 + U — h — I2 + 2 are 
zero since in those cases the respective numerators vanish. Alternatively we can integrate over the hypergeometric 
function directly and write the solution in the following closed form 



A3 AT 



1 



1 



34 ■HihLhisUL 2^ ^2^^ + i)(2l2 + l){2h + 1)(2?4 + 1) 



i?(L + 1,1/2) 



1 



1 

1 1 



B 



1 



1 



E'.+4 A2 

1 1 



B 



B 

C_2 

2 
B 



Ci + 2 



, 1 3^^2 { 



-Ci+2 1 



3 Ci+4 



2 ■'2'^+^>'{^+2' 2 



,l)3^^2({^,^,i + l};{L+^,^};l 



C2 



, 1 3i^2 { 



C2 + 2 1 



3 C2 



2 -'2'^ + l}={^ ■ 2' 2 



1 I , (-7 1 1 ..^ c7 1 -f- 2 ^ 



(70) 



where B[x, y) denotes the beta function and Ci = 1 + Z3 + Z4 + L, C2 = 1 + + /2 + 

The equilateral shape has also received a lot of attention in the literature. As has been described in [T3], for the 
purposes of data analyses, there are two representative forms for the equilateral trispectra. These are given by the 
following shapes for the reduced trispectra 



5;9^|\ki,k2,k3,k4;K) cx K 



c,cquii/, , 1 1 T^N kik2hkiK 
■^r si (^i''*2,k3,k4;K) cx - — — — — 



fcifc2fc3fc4 
(fci + fc2 + fc3 + ki)^ 

2 / 1 



6(fc3 + fc4 + KY 



(fc3 + fc4 + Kf \2{ki+k2+Kf (fci + fc2 + fc3 + fc4)S 

^ 3(fc3 + k4 + K) ^ 1 

(fci + fca + fc3 + kiY (fci + fc2 + fc3 + fc4)^ 



(71) 



(72) 
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where we use the notation ci and si to correspond to [T3]. These shapes are similar in most regions apart from the 
doubly squeezed limit {k^ = ^4 — 0). It has been observed that the first ansatz is factorisable by introducing the 
integral 1/M" = (l/r(n)) f^-^e'^^dt where M = Y^h. As we observe from Figures [g] and [7| it is clear that the 
shapes for the two representative forms are highly correlated. Therefore, for the purposes 01 analysis of the equilateral 
model, it may only be necessary to consider the ci model. 



V. MODE DECOMPOSITION 



Our goal is to represent an arbitrary non-separable reduced primordial trispectrum (zero mode) 
7'^^o{ki,k2,k3,ki; K) or reduced CMB trispectrum Ti''^i^{L) on their respective wavenumber or multipole domains 
using a rapidly convergent mode expansion. We need to achieve this in a separable manner, in order to make 
tractable the five dimensional integrals {'^ dkidk2dk^dkidK) required for trispectrum estimation by breaking them 
down into products of one-dimensional integrals. In particular, this means that we wish to expand an arbitrary 
non-separable primordial (reduced) shape function in the form 

Sr{ki,k2,k3,k4,K) ^ ^ aprsuvqp{ki)qr{k2)qs{k3)qu{k4)rv{K) (73) 

where the qp, are appropriate basis mode functions which are convergent and complete, that is, they span the space 
of all functions on the wavenumber domain. The differing notation, g, r is due to the different ranges of the variables 
- ki £ [0,fcinax] but K S [0, 2fcmax]- In the case of more general probes of non-Gaussianity this is easily extended to 
include the other Legendre modes of equation (pTj) by writing 



Siki,k2,k3,ki,K,ei) ^^Srn{ki,k2,k3,ki,K)P„{cosei) (74) 



where S is the shape function applied to the full Legendre expansion (21 1. The shape function of the nth Legendre 



mode, Sj-n, may be decomposed as in equation (73 1 



We will present one method for finding the basis functions q, r below. We will achieve this objective in stages. 
First, we create examples of one dimensional mode functions which are orthogonal and well-behaved over the 
full wavenumber (or multipole) domain. We then construct five dimensional products of these wavefunctions 
qp{ki)qr{k2)qs{k3)qu{k4)Tv{K) — Qm- This creates a complete basis for all possible reduced trispectra on the given 
domain. By orthonormalising these product basis functions Qm — Ti-m we obtain a rapid and convenient method 
for calculating the expansion coefficients ap^suv (or Oim)- Here we use bounded symmetric polynomials as a concrete 
implementation of this methodology. Of course as outlined in the case of the bispectrum in [T] there are alternatives 
to using the polynomials Qm,T^m but the shortcoming of these alternatives is that either (i) they can lead to over- 
shooting at the domain boundaries or (ii) the choice may compromise separability. However, it is possible that an 
alternative to the polynomial expansion may be desirable to improve the rate of convergence. These should be able 
to conveniently represent functions in a separable form, and should be derived explicitly for the domain. 



A. Domain and weight functions 

In Fourier space, the primordial reduced trispectrum zero mode 7$, 0(^11 ^27 ^3; k^; K) is defined when the wavevec- 
tors ki, k2, K and ka, k4, K close to form triangles subject to ki + k2 -I- K = = ka + k4 — K. Each such triangle 
is uniquely defined by the lengths of the sides ki,k2,K and k3,ki^K. In terms of these wavenumbers, the triangle 
conditions restricts the allowed combinations into a region defined by 

fci < if + /j2 for ki > k2,K, or kz < K + ki for ks > ki, K, or K < ki ka for K > ki, ka (75) 

and 

k3<K + ki for kg > k^, K, or k4 < K + kg for k4 > kg, K, or K < kg + k4 for K > kg, k4. (76) 



Since the wavenumber K is common to both triangles the region is a product of the tetrahedral domains defined by the 
conditions ( 75 1 and ( 76 ) . Considering each region individually we note that they each describe a regular tetrahedron 




Figure 6: Equilateral ci model ( 71 1. The signal peaks at e = 1 towards 7 = 0, a = 0, /3 = 0, i.e. at k\ = k2 = = = K/2. 
For e > 1 the signal similarly peaks for k\ — fe, fea — ki but since the triangles (fci, ^2, K) and (fca, k^, K) are now unequal the 
peak position is less sharp and shifts to smaller values of K. 
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Figure 7: Equilateral si model (72 I. As for the equilateral ci model the signal peaks for equal ki at e, while for e > 1 the peak 
position becomes less sharp. As is clearly observable from the figures the equilateral ci and si models are highly correlated. 
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for ki + k2 + K < 2fcmax (or + k4 + K < 2fcmax)- However, motivated by issues of separability and observation, it 
is more natural to extend the domain out to values given by a maximum wavenumber fcmax- In particular, we have 
ki < fcmax and K < 2fcinax- In each case the allowed region is a hexahedron formed by the intersection of a tetrahedron 
and a rectangular parallelepiped. For brevity we will denote this configuration as a tetrapiped. This region is an 
extension of the tetrapyd referred to in [T] due to the extended range of K and is shown in Figure Is] 




Figure 8: 'Tetrapiped' domain for allowed wavenumbers of the primordial reduced trispectrum T(/ci, k2, ks, /s4; K, 64) imposed 
by the triangle created by (fci, ^2, K). There is a corresponding tetrapiped domain imposed by the triangle created by (fca, ^4, K). 
The region is an extension of the tetrapyd domain described in [T] due to the extended range of K. The same domain is valid 
for allowed multipole values U,L in the case of the reduced CMB trispectrum Ti^i^{L). The shaded area denotes the region 
described in Figure [3] 

In order to integrate functions /(fci, k2, k^, k^, K) over the tetrapiped domains, which we denote V7-, we note the 
presence of K in both regions. We find that the integration is given explicitly by 

![/] = / f{kuk2,k3,ki,K)w{ki,k2,k3,k4,K)dVT (77) 

pl/2 / pt pt+s I'l-t pt+s pi \ 

= fcj^ax \ dt \ ds dx + ds dx + ds dx \ 

X [ dy dz+ dy dz + dy dz) FW 

\Jo Jt-y Jt Jy-t Jl-t Jy-t J 




where uj(k\^k2-,k3,k/^^K) is an appropriate weight function and we have made the transformation i ~ Kjk^^^^s = 
ki/kyne.x,x = fc2/fcmax,y = h/k„ia.x,z = k^/k^^^ with F{s,x,y,z,t) = f{k„-,^^ X {s,x,y,z,t)) and W{s,x,y,z,t) = 
^(/cmax X {s,x,y, z,t)). FoT integrals over the product of two functions / and g we can define the inner product 
{fi9) = 2^1/5] • This inner product essentially defines a Hilbert space of possible shape functions in the domain. The 
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total volume of the domain is given by = fc^ax/^- Initially we will restrict attention to the case of weight uj = 1. 
However, it is important to incorporate a weight function for a variety of reasons which we will discuss later. 

Analysis of the CMB 'extra'-reduced trispectrum, t\]j^^{L), is more straightforward than in the primordial case. 
This is because the CMB trispectrum, being defined on a two sphere, is an explicitly five dimensional quantity and 
therefore is defined completely in terms of the multipoles. We note here that the quantity pJ^^(L) = + 
+ tj^l^(L) + t|^|^(L) is probably a more elegant expression for this analysis since it is more symmetric, whilst 
being defined on the same domain and being subject to the same weighting over the domain. Nonetheless, we proceed 
in this paper with analysis of leaving exploration of this minor issue to an upcoming paper |41| . As for 

the primordial case, we extend the tetrahedral domains to include multipoles out to < /max- The respective 

tetrapiped domains for the extra-reduced trispectrum becomes the discrete {h, hihih, L} satisfying 

h < h + L for li > li,L, + cyclic perms 
h ^ U + L for I3 > I4, L, + cyclic perms 

h+l2 + L = 2ni, l3 + li + L^2n2, ni,n2eN. (78) 
In multipole space, we will be primarily dealing with a summation over all possible {Zi, ^2, ^3, ^4, combinations in 



the correlator C(T, T'). The appropriate weight function in the sum from (23) is then 



uj{h,l2,h, k, L) = hj^i^^h^i^^. (79) 

A straightforward continuum version of this can be deduced by comparison of this 'weight' formula to the bispectrum 
multipole weight function in [T]. Similarly to that analysis, we should eliminate a scaling in this weight such that the 
overall weight becomes very nearly constant. We can do this by using a separable weight function as 

,,(1111 r^ - ^{h,l2,h,k,L) 

' - (2^^ + 1)1/3(2/2 + 1)1/3(2/3 + 1)1/3(2/4 + 1)1/3(2L + 1)2/3 ' ^^"^ 

We note that there is also a freedom to absorb an arbitrary separable function vi into the weight fucntions. If we 
define a new weight uj in the estimator as 

^iii2hUL ^ uji^i^i^ulI {vitVi^vi^vi^VLf , (81) 

then wc must rescale the estimator functions by the factor vi-^vi^vi^vi^vl- The important point is to use both the 
weight UJ and the estimator rescaling throughout the analysis, including the generation of appropriate orthonormal 
mode functions. 



B. Orthogonal polynomials on the domain 

We now construct some concrete realisations of mode functions which are orthogonal on the domain Vj- and which 
have the form required for a separable expansion. First, we will generate one-dimensional orthogonal polynomials 
qp{s),ry(t) for unit weight uj ~ 1. Considering functions qp{s) depending only on the s-coordinate"^ , wc integrate over 
the a;, ?/, z-directions to yield the weight functions uj{s) for s G [0, 1] (for simplicity we take fc,nax = 1): 

5 fi 

uj{s) = s-s^ + — s^ with = / f(s)w(s)ds. (82) 

12 Jo 



We can consider s as corresponding to any of the ki. 
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Therefore, the moments for each power of s become 



1 



1 



n + 2 12(n + 5) n + 4 

5n2 + 54n + 160 
12(n + 2)(n + 4)(n + 5)' 



(83) 



For functions r„(i) (where i corresponds to the K coordinate), we integrate over the s,x,y,z- directions to yield the 
weight functions u>{t) for t G [0,2]: 



^(t) = (^(4-30) forte [0,1] 



forte [1,2], with X[f] = / f(t)t;;(t)dt. 



With this choice of weight the moments of each power of t become 

n2 + 15n + 68 768 x 2" - 744 - 474n - ISln^ - IBrr^ - 

"~ ^ ^ 4(n + 3)(n + 4)(n + 5) 4(n + l)(n + 2)(n + 3)(n + 4)(n + 5) 

Prom these moments we can create orthogonal polynomials using the generating functions. 



(84) 



(85) 
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(87) 



where we choose the normalisation factors J\fi , A/2 such that X[(j'^] = 1 and X[r^] = 1 for all n e N, that is such that 
the qn{s) (or r„(t)) are orthonormal 



(9n,9p) =2:[g„gp] = / qn{s)qpis)dVT = S„p, 
Jvt 

{ry,ru) =I[ryru] = / r„(t)r„(f)dVr = <^^«- 



(88) 
(89) 



The first few orthonormal polynomials on the domain Vr are explicitly 

qo{s) = 

qi{s) = 7.16103(-0.608333 + s) 

g2(s) = 7.76759 - 33.2061s + 29.0098s^ 

g3(s) = -11. 7911 + 93.1318s -194.111s2 + 116.964s^... 



(90) 
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and 

roit) = 

ri{t) ^ 6.06977(-0.7 + t) 

r2{t) = 7.65066 - 24.1493i + 16. 1942i2 

rsit) = -12.2182 + 63.4315t- 91.6438^2 +38.7091^^... (91) 

We note that the g„,r^,'s are only orthogonal in one dimension (e.g. ((J'n(s)r^,(i)) ^ and {qn{s)qm{x)) ^ 5nm in 
general). However, as product functions of t, s, x, y and z they form an independent and well-behaved basis which we 
will use to construct orthonormal five-dimensional eigenfunctions. In practice the g„'s and r^'s remain the primary 
calculation tools, notably in performing separable integrations. In Figure |9] we plot the first few g„ and r„'s. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



60 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 



Figure 9: The orthonormal one-dimensional eigenmodes qn,rv plotted on their respective domains for n,v — 0,1,2,3. The 
qn,rv plotted are calculated for unit weight = 1 on the domain Vt- The shape of these eigenmodes alters for different choices 
of the weighting. 



24 



Now we turn to the polynomials g(s) and f(t) which are orthonormal on the multipole domain. Using the weight 
function uj = 1 we find the same polynomials as above. For the scaled weight function Wg the polynomials will, of course, 
differ. While using either polynomial sets would suffice as independent basis functions on the multipole domain, the 
use of correctly weighted functions lead to improvements in the immediate orthogonality of the derived five dimensional 
polynomial sets. For definiteness we take {t,s,x,y,z) x Imax = {L, 11,12,13,14). The generating function is obtained 
as above but now using the moments w„ = I[s"] — J oj{t, s, x, y, z)s"c?Vr and a)„ = = J uj{t, s, x, y, z)t"dVr- 



C. Five-dimensional basis functions 



We can represent arbitrary (reduced) trispectra Legendre modes on the domain Vj- using a suitable set of in- 
dependent basis functions formed from products qp{s)qr{x)qs{y)qu{z)r-u{t). (Here again we take t = -fC/femax, s — 
^i/fcmax, X = fc2/fcmax, V = ^s/^max, z = k^/kma.^ OT t = i/^max, etc.) Wc denotc the 5D basis function as 

Qjn{t,s,x,y,z) ^ qp[s)qr{x)qs{v)qu{z)ry{t). (92) 

We can order these products linearly with a single index m in a similar manner to that described in 1 for the 
bispectrum. 

While the Qm's by construction are an independent set of five-dimensional functions on the domain V7-, they are 
not in general orthogonal. To construct an orthonormal set TZm from the Qm we perform an iterative Gram-Schmidt 
orthogonalisation process such that 

In particular, we form the Gram matrix F = {{QnQm}) which needs to be factorised as F = A-'" A where A = {{QnTlm)) 
is triangular. This process is described in more detail in [T]. 



D. Mode decomposition of the trispectrum 

Having formed the orthonormal basis {TZm} we consider an arbitrary primordial reduced trispectrum (zero mode) 
Tfl{ki, ^2, fcs, ^4; K) described by the shape function Sj- and decompose it as follows 

Sr{ki,k2, ks, ki, if) = ^ a'^Tlm{t, s, x, y, z) (94) 

in 

where the expansion coefficients are given by 

= {n,n,ST)= ! n^nSrUjdVT (95) 

Jvt 



c3 



and fcmax(^, s, X, y, z) = (K,k\,k2,kj,,k^ on the domain V7- defined in (75) and (76). In practice, we must always 
work with partial sums up to a given N = rtmax with 



=Y^a^Tlr,^{t,s,x,y,z), 5r = lim 5^. (96) 

^ ^ 7-00 

Given the complete orthonormal basis 7?,^, Parseval's theorem for the integrated product of two functions implies 

{St, S!r) - / SrS'rUjdVr = lim ^ a^a^ (97) 

which for the square of a mode Sj- yields the sum of the squares of the expansion coefficients, I[S'^] = J2m '^m ^■ 
In order to accomplish the goal of a general separable expansion, we must transform backwards from the orthogonal 
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sum TZm into an expansion over the separable product functions Q„j through 

N 

=Y.a^Q^{t,s,x,y,z). (98) 

The can be obtained from the via 

N 

- E(^^)™f"? (99) 

where the transformation matrix A„p was defined above. Using the inverse relation — X]^o('^"^)mp'^^ f^'^^ 
that the matrices F and A are related by 

N 

(7-^)„,-E(^")„.^^.- (100) 

r 

This implies that 

N N N 

(^f > 5f ) = E "™ ' = E E "^.7™pa?. (101) 

m—Q rn—0 p—0 

As we have already noted the separable Qm expansion is most useful for practical calculations. However, its coefhcients 
must be constructed from the orthonormal TZm- 

We can expand the CMB extra-reduced trispectrum i;|^|^(L) at late times using the same polynomials. However, as 
noted previously the CMB trispectrum is an explicitly five dimensional quantity and, as such we do not require the 
extra step of finding the zero mode of the Legendre series expansion. In particular, the appropriate expansion is of 
the form = J^m'^mT^mit, s, x,y, z){= T„n^niQr,i{t,s,x,y,z)). 

VI. MEASURES OF Tnl 
A. Primordial estimator 

We have obtained related mode expansions for a general primordial shape function, one with the orthonormal basis 
TZm and the other with the separable basis functions Qm- Substitution of the (reduced) separable form into the 



expression for the 'extra'-reduced trispectrum ( 23 ) offers an efficient route to its direct calculation through 

t\l\l{L) = NAl 0^ J rldrirldr2dkidk2dks,dkidKK (102) 
xE"™2m(fci,fc2,fc3,fc4,i^)jL(^ri)jL(ifr2)[j,,(/ciri)Ai,(fci)][j;,(fc2ri)A,,(A:2)] 

m 

= 7VA|E«^ / rldr,rldr2Q'^''''^'^{r^,r2) 

where 

Qji'^'^'-^(ri,r2) = ql,'{ri)qHri)qHr2)qlt{r2)r^{ri,r2) (103) 
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with 



i{r) 



dkqp{k)Ai{k)ji{kr), 
dKKr,{K)jL{Kr^)jL{Kr2) 



Next, we note from (11) that 



LM 



U I9. L 



M / '1 '2 



13 14 



L L 



mi 7712 —A/ / V 7713 7714 M J '3^4 ' 



(104) 



(105) 



Therefore, using these formulae in the estimator (28) (we omit the normaUsation factor here and return to it later 
in the section) we find 



lirrii 

£ = ^^12A^A 



Ci^Ci^Ci^Ci^ 



lim: LM 



dniYi^^^{hi)Yi^rn2ini)Ylj^,j{hi) j dn2Yi^^.-^{n2)Yi^rnAf>-2)YLM{n2) 

(^^imi ^l2m2 ^l^m^ ^141714 )c 



(106) 



Now using the notation £ — £tot — ^uc where tot refers to ai-imiO'i2m2'^i3m30'iimi and uc refers to 
(a; )^, we find 

£uc = 12N a^^ J dhidn2 J dridr2rlrlN^{ni,n2,ri,r2)^Mp^{hi,hi,ri,ri)M^^{h2,^ 

+M^^{ni,n2,n,r2)M^^{hi,h2,ri,r2) + M^^{hi,h2,ri,r2)M^^{n^ I (107) 



and 



£tot = 12NA%'^a^-^J dhidn2 J dridr2rjrlMp{ni,ri)Mr{hi,ri)Ms{n2,r2)Muin2,r2)Ny{hi,h^ 



(108) 



where 



Mp^{ni,n2,ri,r2) = 2^ 



1 1 mi 



Mp(ni,ri) = J2 



1 1 mi 



Ny{hi,n2,ri,r2) = ^Yl^.j{hi)YLM{n2)ry {ri,r2). 



(109) 



LM 



We can summarise these results (substituting back in Nt) as 

\2NA^ > f f 

£ = — — — ^^"m / dnidn2 / dridr2rlrlM^{hi,h2,ri,r2) 



Nt 
Nt 



e 



(110) 
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with 

(3^ = 12 J dhidh2 J dridr2rlrlM^{ni,n2,ri,r2) (111) 

and the form of Ai^^ inferred from the above equations. The estimator has been reduced entirely to tractable integrals 
and sums which can be performed relatively quickly. 

We can estimate the computational time needed to find this estimator. The multipole summation needed for each 
basis is 0(?inax) (since the sum over the m's can be precomputed) . The integral ^ J d,^f>- is an C(^max) calculation, 
while the line of sight integral ~ J dr is conservatively estimated as an O(IOO) operation. Therefore, in total the 
estimated number of operations is 0(10000) x 0(1^^^.). 

In the case that the primordial trispectrum is independent of the diagonal K, the estimated number of operations 
reduces to O(IOO) x 0(/i^ax) ^ outlined in Appendix D. 



B. CMB estimator 



In the case of a precomputed CMB trispectrum or a late times source of non-Gaussianity in the CMB, such as 
gravitational lensing or active models such as cosmic strings, we wish to find a late-time CMB estimator. For the 
late-time analysis we wish to expand the estimator functions using the separable Qm(^i, '2, ^3, ^4, i) mode functions 
created out of the qp{l) and fy{L) polynomials (Note that we denote the multipole modes with a bar to distinguish 
from the primordial equivalents, and also that we have no need for a subscript for the zeroth Legendre mode since 
the CMB trispectrum is an explicitly five dimensional quantity as described earlier) . In order to effectively expand in 
mode functions modulated by the C;'s we choose to decompose the estimator functions directly as 

^;tfff^4!:(i)-E=sa,. (112) 

where the separable vi incorporates the freedom to make the weight function uj even more scale invariant. The 



estimator expansion with Ci in ( 106 1 is appropriate for primordial models, but it is expected that flatter choices will 



be more suitable for late-time anisotropy, such as from cosmic strings. 



Substituting this mode expansion into the estimator ( 28 ) (where again we omit the normalisation factor Nt and 
return to it later in the section) we find, 

^ = ^'2^^^a^9p{h)qr{h)qs{l3)qu{h)rv{L) / d^niFiimi (ni)!!^™, (fii)Fj^V(^i) 



lirrii LM n 



X / d'h2Y,,^,in2)Yi,^An2)YLMin2) ' (113) 
After some algebra we find 

5*°' - 12ArA|^cF2 f d^hid^h2Mpini)Mr{ni)Msih2)Muin2)Nvin^^ (114) 
f"'^ = 127VA|^a2 f d^hid^n2 

n •' 

"lie "lie "LLC UC UC "UC \ 

Mp^{hi,hi)M^^{h2, ^2) + Mp^{ni,n2)M^^{ni,n2) + Mp^{ni,n2)M^^{ni,h2)j iV^,(ni, ^2), 



28 



where 



lirrii 

mi 

ih)Uli) 



Again we can summarise these resuhs (substituting back in Nt) as 

/VA; 



with 



= 12 / d^nid'h2M'^ini,n2), (117) 



and the form of Ai^ {ni,n2) can be deduced from the equations for and S^'^. 

Since there are no Une of sight integrals (~ J dr) for this estimator the number of operations required in this case is 
^Cmax) suggesting that the late-time estimator is much more computationally efficient than the primordial version. 

Similarly to the primordial case, there is a reduction in complexity to ©(/^ax) in the case that the late-time 
extra-reduced trispectrum is independent of the diagonal L. This is explained further in Appendix D. 



C. TjvL Estimator 

As with shortcomings of normalising the quantity Jnl of the bispectrum that was addressed in [1] the current 
method |14j of normalising the level of non-Gaussianity due to the trispectrum, tj^^ poses problems. In particular, 
the level of non-Gaussianity is found by normalising the shape function against a central point. More specifically 



we can identify this method as setting Sf{k,k,k,k,k) = 1 and identifying the normalisation N of equation (37 1 as 
(50/27)tArL. In the case of the local model this gives 

t'^t - 2.16/^i = l.SrjvL , t'°f = l.OSgNL, (118) 

where we note again that the relationship between tjvl and /^vl is only strictly true for single field inflation. This 
approach assumes scale invariance and therefore will produce inconsistent results between models peaking or dipping 
at this central point. Also, this definition is not well-defined for models which are not scale- invariant, such as feature 
models, and it is simply not applicable to non-Gaussian signals created at late times, such as those induced by 
cosmic strings or secondary anisotropies. An alternative measure of the non-Gaussianity is given by comparison of 
the primordial trispectrum to the local primordial trispectrum but this approach is not well-defined and is essentially 
an order of magnitude estimation jl3j . 

Therefore, we propose a universally defined trispectrum non-Gaussianity parameter T/vl which (i) is a measure 
of the total observational signal expected for the trispectrum of the model in question and (ii) is normalised for 
direct comparison with the canonical local model (with g^^ — 0). We define Tjvl from an adapted version of the 



estimator ( 28 1 with 



1 /n n n n \ ( n'^^^ ^obs „obs ^obs \ 
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where N is the appropriate normahsation factor for the given model, 



and -/VtiocA is the normahsation for the focal model with tj^l = 1, 9nl = 0. 



The Tatl estimator will recover ttvl for the local model with gjvL = 0, while it gives {Ntiocb/ NT\ocA)gNL for the local 
model with tj^il = where -/VtiocB is the normalisation for the local model with rjv^ = 0, gjqL = 1- This coefficient 
is dependent on Zmax but is a number of order unity. 

Results for primordial models should not depend strongly on the multipole cut-off /max- However, diffusion due to 
Silk damping in the transfer functions ensures that the primordial signal is exponentially suppressed for I > 2000. 
Therefore, an appropriate choice for a canonical cut-off is /max — 2000. Late-time anisotropies, such as cosmic strings, 
do not generically fall off exponentially for / > 2000 but, nonetheless, in the domain I < 2000 we can make a meaningful 
comparison to the local tnl = ^,gNL = model. Alternative measures must be proposed beyond this domain. As 
indicated in Appendix A the normalisation factor is a computationally intensive calculation. Instead we use the 
approximation Nt ~ A'^tiocaIAV/A^TIoca) where the subscript T instead of T refers to using the reduced trispectrum 
instead of the full trispectrum in the above calculations. With these approximations we need only accurately calculate 
the full normalisation factor in the case of the model with tnl — ^,gNL — 0. Regardless of the accuracy we propose 
that given the vastly increased speed of the calculation we adopt this latter convention and define T^l^ i.e. 

-,obs ^obs ^obs ^obs 



/'2^yoc(rNL = l,gNL=0) 



~\f I n n n n \ I n^^^ n^^^ /^'-'^^ 1 

rp JVtIocA \^ \"iinil"'2m2"i3m3"'4ni4/c V"iimi"i2m2"i3m3"/4in4jc rinn\ 

J-NL- =2 n n n n • ^^^^> 



The relation between T^l and T^l as well as the accuracy of the above approximation for Nt - which is only a 
conjecture at present- will be explored further in an upcoming paper. However we note here that the Tat^ estimator 
also recovers t^^l in the case of the local model with g^L — 0. 

If the CMB trispectrum is not known precisely for the primordial model under study, then we can make an estimate 



for the normalisation factor A't- in (122) using the shape function for the reduced trispectrum S-j-{ki,k2,k3,k4, K). 



One can obtain a fairly accurate approximation to the relative normalisations in (122) from 



= F{St,St)= j dVkS^r{kiMM,ki,K)u{kuk2MM.K) (123) 



where the appropriate weight function was found in ( 47 1 and the domain Vk is given by ( 75 1 , ( 76 1 . Using A^ / A^riocA ~ 
A^/AiocA to approximate A^t- we can make a fairly accurate estimate of the level of non-Gaussianity or can renormalise 
Tnl constraints for different models into compatible constraints in a similar manner to the analysis of /jvl constraints 
in the case of the bispectrum in [T]. 



VII. RECOVERING THE TRISPECTRUM 
Recovering the primordial trispectrum 



The form of the estimator in ( 106 1 suggests that further information may be extracted from the observed trispectrum 
beyond the tjvl for one specific theoretical model. This is because, through the coefficients /3® we have obtained 
some sort of mode decomposition of the trispectrum of the observational map. Consider the expectation value of 
obtained from an ensemble of maps generated for a particular theoretical model with shape function Sj- = a® Qm- 
Since the shape function is in terms of the zeroth mode of the Legendrc expansion of the primordial trispectrum, we 
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can only hope to recover information about this mode via recovery of the shape function*. Using the expression 

(^J dniYi^rni{ni)yi2m2{ni)^LM{fi-i)j i^j dh2Yi^rn3{n2)yum4{n2)YLM{n2) 



ws.) = 12 EE 

Umi LM 



I'' 



dri4dr2QT''''^{ri,r2) 



' 

as well as the identity for the Wigner 6j symbol in Appendix A we find, after some algebra 

li,L m' V 



(124) 



L+L' I h h L 



dridr2rlriv'J''''''''{n,r2) 



(125) 



where P„ 



The quantity Tmm' represents a matrix with positions 



labelled by m, m' and can be inferred readily by the above equation. Inverting the relationship we can recover the 



aS, via 



c^S = E(r"')w(/3^,). 



(126) 



Therefore, if the decomposition coefficients are found with adequate significance, we can reconstruct the shape function 
through the expansion 



(127) 



This reconstruction will be sufficient to uniquely define the planar case {64 = 0) (as well as the general CMB case in 
the next section) . However, as already discussed the shape function only gives information about the zeroth Legendre 
mode of the primordial trispectrum. Therefore, recovery of the full primordial trispectrum is compromised by this 
degeneracy. However, as we have discussed in Section H, this degeneracy may be broken by using other probes of 
non-Gaussianity such as galaxy surveys and 21 cm observations. We note also that the calculation of the matrix F is 
computationally intensive due to the presence of the Wigner 6j symbols. Nonetheless we include the discussion here 
for completeness. 



^ This is to be somewhat expected since the primordial trispectrum is a six dimensional quantity whereas the CMB trispectrum is explicitly 
five dimensional. 
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Recovering the CMB trispectrum 

The recovery of the CMB bispectrum from a given observational map is more straightforward (as for the bispectrum) . 
In a similar fashion to the calculation in the case of the primordial trispectrum we find that 



(a5 



lirrii LAI 



X / r^drir2dr2Q.ri \ri,r2) 



Vh^h^is w/4 vl Ci^Ci^Ci^Ci^ 



U /.s L' ' P 



Vr 



(128) 



where q'"'''''"^ = gp(;i)g^(/2)g,(/3)9„(^4)r„(L) and T'^i'^'^'"^ = gW2/3UL^_gWii3i4L + gW2Ui3L + gWi;4i3L^ Inverting 
the matrix r„p we find 



«,? = Er»p(^p)- (129) 
— Q 

Therefore, if we can measure the coefficients /3p with significance from a particular experiment, we can reconstruct 

^ Ci^Ci^Cifii^ 



the trispectrum map using (112) 



(130) 



The calculation of the matrix F remains computationally intensive but it is tractable. We can, in principle, extract 
the full CMB trispectrum which, together with the extracted CMB bispectrum [1 , will prove to be a key test of the 
Gaussianity of the Universe. 



VIII. MAP MAKING 

In this section we derive an algorithm for creating a non-Gaussian map with given trispectrum, developing methods 
presented for the bispectrum in ref. |42j and generalised in ref.[T]. This algorithm is valid in the limit of weak 
non-Gaussianity. 

We define the function 



G 



(131) 



LiTTli 



where is the Gaussian part of the CMB multipoles, generated using the angular power spectrum Ci, while 
Ti^rnii2m2i3m3Um4 IS the givcu trispcctrum of the theoretical model for which simulations are required. 
Setting 



*G *G 



1 d 



+ ^ E ^"2(3 



C12 Ci, 4 5a;,„ 



T2[C~'d 



*G *G 
2/3 hm2 "'13^3 

rnm2m3 
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^l2m2 '^l3m3 ^Ijmi 
I2 



(132) 
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we recover the bispectrum from {a.[^^^a'i^^^a'i^^,^) (as described in [35]). Next, we calculate the four point correlator 
of the a',^ 's and find 



a^'l a1^^ al^ 



II I I I \ _ iG „G „G „G \ I /JL T "-ifcmt, "-l^m^ "-Idrng G „G „G , 

ij-mj ''^ 

+perniutations, (133) 

where j G {b, c, d). The first term clearly gives the unconnected component of the four-point correlator. We note that 
the contribution from the bispectrum is zero since the correlator over an odd number of vanishes. Evaluating the 
correlators in the second term on the right hand side and adding up the different permutations we find 

I I I I I \ I Gr G (jr C? \ ' I ' / \ 

\^limx'^l2m2'^hm3'^Um4) ^ ("^hnii '^i2m2 ^'ams '^Uni4 / Ti-^mil2m2l3m3Urn4,- (134) 

This verifies the validity of the use of aj^ to make maps including the Gaussian, bispectrum and trispectrum contri- 
butions to the model under study. 

We observe, using (11 1, that we may write T2[a^] in the form 

T2[a^] ~ 2 '^I™l'2ni2'3ni3'4m4^hmi'^i^)n2'^i3m3'^'4ni4- (135) 

liirii 

Using this formula we may also rewrite the trispectrum contribution to a[^, which we denote a^''^ as 



NG,T ^ _\^T "i2m2 "iama "i4m4 oR\ 

— „ . / , lml2rn2l3m3Um4, f~t K^'-''^) 



nr. 



2^ J ^ t7;tt2 ''t2t3 //t3t4 "t4 



*G *G *G 



_ ^\^(n- in- -i_ T \ "i2™2 "'hm3 "'Um 4 

Q / , \ I I'ml2m2l3m3l4m4 ~r i I2ra2lml3m3l4m4 ~r i I2m2l3m3lml4m4 ~r i I2ni2l3m3l4m4lm } 

'^'3 



Using the formulae for the extra-reduced trispectrum ( 23 ) and the Gaunt integral ( 14 ) , we note that 

Tl,m4l2m2hm3l4m4 = ^^^"1 (ni)^i2m2 (nO^^iV ("l )^ (^J dnn^Yi^„,^{h2)Y^^^{h2)YLM {n2)^ t\l\l{L) . 

(137) 

As an aside, we note that if i||^|^(L) is independent of the diagonal L then we can use equation (171 ) to write 

Ti^mii2m2i3m3i4m4 = / df^n^limi (n)l/2m2 (n)y;3m3 (n)!^^™^ (n)i^|^ , (138) 



where we drop the label L from the extra-reduced trispectrum. This special class of trispectra is explored further in 
Appendix D. 

Denoting the bispectrum contribution to aj^ as a^^'^ we have verified the following prescription for forming maps 
including the bispectrum and trispectrum contributions 

/ G , r ~NG,B , ~NG,T /i on\ 

O-lni^ + JNLai^ +TNLai^ (139) 

where we have written af^'^ — Inlo,^'^ and a^'^ — tatlo^'"^ to make the size of the respective non-Gaussian 
components more explicit. 

Since the computation of the reduced trispectrum is more efficient using the late time expression (due to the 
absence of the line of sight integrals) we write out the formula for a^''^ using the late-time mode decomposition. It 
is straightforward to find the equivalent formula using the primordial expression. Later in the section we present a 
particular application using the primordial local model of this formalism. 

The latc-timc mode decomposition of the extra-reduced trispectrum i|^^(L) - as detailed at the end of section V - 
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may be written as 

Using these expressions we have 



(140) 



li'trii 



where 



„*G „*G „*G 

"'i2"i2 iams "'i4m4 
C12 C'/4 



(141) 



A^r(ni) 



r,,„,(ni)a,*' 



E'l2'"2\---L/"J2m2- /; \ 
77 9rl'2j, 

i2»n2 



C;2 



7V„(ni,n2) = 7V«(n2,ni) = ^y£A^(ni)yi,Af(n2)r^(i). 



(142) 



Evaluating, in a similar way, the other terms in equation (136) wc find 

NG,T 1 



n 



(143) 



As emphasised in [T the condition that the map has the power spectrum C; specified in the imput will only be satisfied 
if the power spectrum of the non-Gaussian components is small. Therefore, one has to ascertain that spuriously 
large C/^*^ contributions do not affect the overall power spectrum significantly. We will discuss the implementation 
of the algorithm presented here in an upcoming paper 41,. 



Application to the Local Model 

The reduced trispectrum for the local model, as shown in Section IV, is made up of two terms which we denoted 
locA and locB. As shown in [33] - and can be deduced from Section IV - the extra reduced trispectra may be written 
as 



where 



t\\\:{Lr^ = 9NL I drr'pi,{r)Pi,{r){piJrMr)+pi,{r)fii^^^^^ 



,locE 



FL{ri,r2) = ^ I K^dKP^{K)jL{Kn)jL{Kr2), 
ai{r) = Mi(^) = ^ / k^dkAi{k)ji{kr), 



A(r) = - / k'dkP^{k)Ai{k)jiikr). 



(144) 
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Using these formulae, and exploiting that the locB trispectrum is independent of the diagonal L with equation (1381, 
we find 

Q*G Q*G 2^ /■ r 

xMain2,r2)Mf3in2,r2), (147) 
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where 
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We similarly find the other terms in ( 136 ) to get 



(»tof'^)iocA = ^"^^i j dridr2rlrl ai{ri) J dni,^dnn^Yi„,{hi)M,3{ni,ri) + (3i{ri) j dni,^dni^^Yi^{ni)Ma{ni,ri) 

(150) 
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(151) 



In the case of the bispectrum, direct implementation of the explicitly separable local shape results in spuriously 
large contributions. However, it was found that using the eigenmode expansion in ref. ^ was much more 

robust circumventing such effects because of the bounded nature of the polynomial eigenmodes. This improvement 
is expected to occur for the trispectrum. An alternative method is to regularise the expressions given here by 
eliminating pathological terms, while leaving the final trispectrum of the map unchanged. For arbitrary separable 
trispectra (unlike the eigenmode expansion), convergence must be achieved by hand on a case- by-case basis. 



IX. CONCLUSIONS 



We have described in this paper two comprehensive pipelines for the analysis of general primordial or CMB trispec- 
tra. The methods are based on mode expansions, exploiting a complete orthonormal eigenmode basis to efficiently 
decompose arbitrary trispectra into a separable polynomial expansion. These separable mode expansions allow for a 
reduction of the computational overhead to tractable levels, regardless of whether the reduced trispectrum is being 
computed at Planck resolution, or we are directly finding an estimator for the size of the trispectrum from a real data 
set. A shape decomposition has been described allowing for a visualisation of a scale invariant reduced trispectrum 
on particular slices. 

We have presented a correlator for comparing trispectra. We have also defined a correlator for comparing the shape 
functions that is expected to closely approximate the former. However, the main purpose of this paper was to present 
a detailed theoretical framework for finding an estimator for the size of the trispectrum using separable eigenmode 
expansions. Using this efficient method for finding an estimator for the trispectrum we have defined a universal 
measure Tjvl which will allow for consistent comparison between theoretical models. This measure can be calculated 
for both primordial models and late-time models, e.g. due to active models such as cosmic strings. The completeness 
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of the orthogonal eigenmodes should allow for a reconstruction of the full CMB trispectrum from the data, assuming 
the presence of a sufficiently significant non-Gaussian signal. We have also detailed an algorithm for producing 
non-Gaussian simulations with a given power spectrum, bispectrum and trispectrum. The implementation of these 
methods will be discussed in a future publication |41j . Clearly, the full implementation of the primordial and late-time 
pipelines represents a significant challenge. However the generality and robustness of the methodology described here 
indicates that there is an intriguing possibility of exploring and constraining a wide class of non-Gaussian models 
using the trispectrum. 
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APPENDIX A: NORMALISATION FACTOR 



Clearly the appropriate normalisation factor for the estimator (281 is of the form 



liirii 



Cl^Cl^Cl^Cl^ 



In order to relate this to Nt expressed in ( 29 ) we use equation ^ to expand Nx' in the form 
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we find that 



{2L + l)Ci^Ci^Ci^Ci^ 



Nt 



This verifies the use of equation ( 29 ) to normalise the estimator ( 28 ) 



We can expand Nt in terms of the reduced trispectrum using 

\ ^ Tlirnil2rn2hnT.3limiTlimil2m2l3m3l4m4 

Then with the identity for the Wigner 6j symbol (see Appendix in [34]) 
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the identities ( 152 ) and relations for P in (|8| and ^ we find that 



(155) 



Due to the presence of the 6j symbols the calculation of Nt is computationally very expensive in general. 

APPENDIX B: OPTIMAL ESTIMATOR 

When non-Gaussianity is weak we can exploit the multivariate Edgeworth expansion [3 7) around the Gaussian 
probability distribution function (PDF), P*^(a), i.e. 



P{a) = 



1 ^ ^ {^linii ^12^12 ^Isriis I 
liirii 



d d d 

ma 



d d d d 

^Tfl4^ / C r\ A O 
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P'=(a), 
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where the Gaussian PDF is given by 



P«(a) 



(27r)^/2|C|i/2 



(157) 



with Cim,i'm' = {aimai'm') and N the number of I and m. Maximising over the three point correlator results in the 
optimal bispectrum estimator. Here we will ignore this term (setting it to zero for convenience) and concentrate on 
the four-point correlator. We find 



P(a) 



1 + (fliimifli2m2'^i3m3'^'4m4)c I (C* "'"a)iimi (C ^0,)l2rn2{C' ^0,)l3m3{C ^(l)l4^ 



P<=(a), (158) 



where (C ^ajim = X^z'm' ^I'm' im^i'm'- Parametrising the size of the trispectrum by £ we wish to maximise the PDF 
with respect to this. We assume that {cii^mx^i2m2^i^m^^i4^m^^c CK S sucli that the second term is proportional to S. 
Maximising the PDF means that we wish to set (dP/d£) = 0, such that the Taylor expansion about ^ = reads 
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Since 

d^P 
d£^ 



141714 ,l'^m'^ 



li-rrii 



we find that the estimator is maximised by setting (with appropriate choice of proportionality constant) 
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where 
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APPENDIX C: NON-OPTIMAL ESTIMATORS - SKEWNESS AND KURTOSIS 



The deviation from non-Gaussianity may be measured in a non-optimal way by estimating the departure of the 
one-point PDF from Gaussian behaviour. This deviation may be measured in terms of the skewness and kurtosis. 
The skewness is given by 



91 = 



3/2 ' 



while the kurtosis is given by (|35|. Using Q we evaluate the variance as 
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The three-point temperature correlator is similarly given by 

3 \ 



.y;„,(n)rz™(n) = i-^(a 



An 



/AT 



) = ^{ahnnai^,n2ai3,n3) J ^^iimi (n)^/2m2 (n)>";3m3 (n) 

/ lit'lli 



where Grnim2m3 is t-iis Gaunt integral and bi^i^i^ is the reduced bispectrum. The Gaunt integral is given by 
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Using equation ( 152 1 we can simplify this expression to get 
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Next in order to evaluate the kurtosis we calculate the quantity 
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where the final three terms are clearly equivalent. Using 

(a/ 37713*^/47714/ "I V'^/iTTli '^/37?l3 

+ {ai iTTli 0-/47714 ), (167) 

it can be shown after some algebra that 
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From equation ( 11 ) this may be written in terms of the reduced trispectrum as 
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Next, noting that the product of two spherical harmonics can be written as 
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and equation ( 14 ) , we find 
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Finally, using equation ( 10 1 for the reduced trispectrum, the orthogonality relation between Wigner 3j functions as 

(172) 



expressed in (18 1 and the 'extra'-reduced trispectrum (|23|, we find 
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In summary, the skewness and kurtosis are given respectively by 
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APPENDIX D: SPECIAL CASE OF TRISPECTRUM INDEPENDENT OF DIAGONAL 

Suppose that the primordial reduced trispectrum is independent of the diagonal K. In particular, we write 



T<s,fi{ki, k2, fcs, k^, K) = 7$,o(fci) ^2? ^3) ^4)- In that case the 'extra'-reduced trispectrum (see (22) and (23l) becomes 



t\±{L) 
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Next, using equation (61) we find 



dKK'jL{Kri)]L{Kr2) = ^5(^2 -ri). 
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This implies that 



( fci fc2 /ca A:4 ) ^ dfci (ifc2 (ifcs (ifc4 r 2 dr 1 



i.e. we only have one line of sight integral. This expression also shows that, if the primordial trispectrum is inde- 
pendent of the diagonal K, then is independent of L. We can exploit this property in our estimators. From 



equations (10), (23) and equation (14) for the Gaunt integral (which we denote here in the form Gm^m^ms) we have 
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If the extra-reduced trispectrum is independent of L we can use equation (|171|) to write 
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where we drop the label L from the extra-reduced trispectrum. The extra reduced trispectrum now has the following 
mode expansion 



where now we have 
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The mode decomposition is similar to that described in section V with r„ — constant. In this case we use the following 
primordial decomposition 



(A:ifc2fc3fc4)^7^^o(fci, ^2, ^3, ^4) = ^ a^iqp{ki)qr{k2)qs{h)qt{k4)- 



(183) 



Using this in the expression for the general estimator (32), which can be re-expressed as 
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we find 



12iVA 



Nt 



dh dr P Mp{h,r)Mr{h,r)Ms{h,r)Mt{h,r) - (Mp^{h,r)Ms{h,r)Mt{h,r) + 5peTins 



- [M^^{n, r)M^[{n, r) + 2 perms 



(186) 



40 



where 

Mp{n,r) = 2^ ^r^9pW, 

Im 

M;^{h,r) = E '^'*"^^^gy""^"^ gp(Og[(r-)- (187) 

We again can estimate the computational time needed to find this estimator. Since we now have only one line 
of sight integral and one integral over the sky ~ J dn we use the prescription outline in Section VI to estimate the 
complexity conservatively as 0(100) x O(^max)- 

The implementation in the case of the late-time estimator for which the extra-reduced trispectrum is independent 
of L can be found similarly. Since this estimator no longer requires a line of sight integral the complexity of the 
calculation can be estimated as O(^max)- 
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